Stable Numerical Integration of an Epitaxial Growth Model with Slope Selection 
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We consider a continuum phase field model for molecular beam epitaxial growth with slope selec- 
tion, with the goal of determining stable numerical time integration methods for the dynamics. We 
parametrize a class of semi-implicit methods that are linear in the updated field, which allows for 
efficient implementation with fast Fourier transforms. We perform unconditional von Neumann sta- 
bility analysis to identify the region of stability in parameter space, and then test these predictions 
numerically for gradient stability. We find strong agreement between the approaches. 
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I. INTRODUCTION 

In molecular beam epitaxial (MBE) growth the 
Ehrlich-Schwoebel effect [l|, Hj can destabilize a flat in- 
terface and lead to the formation of pyramids or mounds 
. These surface features then coarsen, with their height 
and spatial extent growing as powers of time. Theoretical 
studies of MBE coarsening typically employ continuum 
models, justified by appeal to the large distance and slow 
time scales involved. The resulting field equations of mo- 
tion are nonlinear, and to make progress they must be 
integrated numerically, a process which, unfortunately, 
is hampered by numerical instabilities. As such, much 
recent effort has been devoted to finding stable integra- 
tion methods. In this work, we derive a class of stable 
numerical integration methods that are particularly effi- 
cient and simple to implement because the updated field 
can be obtained via the fast Fourier transform (EFT). 

The model we consider employs a height field h(x, y, t) 
that is a continuous function of space and time, and 
which obeys the equation of motion 
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applicable for MBE growth with isotropic slope selection. 
The motivation for this and related models is discussed 
below. With these dynamics, equilibrated regions of uni- 
form gradient and unit slope form. Domains with differ- 
ent slope orientations meet at edges of constant width, 
and as the system evolves the edges are healed out, re- 
sulting in the growth of the characteristic domain size. 
Eor this particular model it has been found from theoreti- 
cal analysis [1, Q , simulations d, 0| , and rigorous bounds 
that the characteristic domain size L(t) grows with 
time as L ~ t^/^. 

Numerical simulations of coarsening are useful for test- 
ing scaling and the predicted growth laws and for measur- 
ing properties of the scaling state, such as correlations, 
growth law amplitudes, autocorrelation functions, and 
more (see Q for a coarsening review). But these simu- 
lations face several restrictions. To reach the asymptotic 



scaling regime, it is necessary to evolve until L{t) ^ w, 
where w is the width of the edges. But the lattice size Aa; 
must be sufficiently smaller than the edge width in order 
to resolve the edge shape and corresponding line tension. 
Finally, the system size Lgys must be large enough that 
domains can grow into the scaling regime before finite 
size effects appear. To satisfy this string of conditions, 
Aa; w ^ L{t) Lsys, requires lattices of very large 
linear size Lgyg/Aa;, evolved to late times. 

For this reason, it is desirable to use integration 
schemes that are accwracy-limited rather than stability- 
limited. Euler integration of Eq. ([1]) is only stable for 
time steps At smaller than a threshold determined by 
the lattice spacing. In contrast, an unconditionally sta- 
ble method, i.e., one with no conditions on At, would 
allow a time step determined by the natural time scale of 
the dynamics, which turns out to be considerably more 
efficient. Accuracy considerations require the typical dis- 
tance traveled by an edge within one step to be held fixed 
and since the characteristic edge velocity scales as 
I'odgo dL/dt ~ t~^/^, this allows a growing time step 
At t^/^. Using dt/dn ~ At, where n is the number 
of integration steps, it follows that unconditionally sta- 
ble methods allow accurate evolution with t ~ n^, rather 
than the stability-limited t ~ n. 

Eyre provided a general approach for generating un- 
conditionally stable semi-implicit integration methods, 
based on a splitting into expansive and contractive terms 
Wang, Wang, and Wise used this approach for 
Eq. Il]), as well as an MBE model without slope se- 
lection and this approach has now been extended 
to a second-order in time method and other de- 
velopments These schemes are gradient stable, 
which means they preserve the energy-decreasing prop- 
erty of the continuous-time equation. However, these 
Eyre-based schemes have the drawback that usually a 
nonlinear term must be treated implicitly, requiring an 
iterative method to find the updated field, and in the 
worst case no guarantee of convergence or a unique solu- 
tion. An alternate approach is to restrict consideration 
to steps with linear implicit terms that can be solved di- 
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rectly by FFT, determine the range of step parameters 
that satisfy unconditional von Neumann (UvN) stability, 
and then test these parameters numerically for gradient 
stability. This approach yielded stable, direct steps for 
the Cahn-Hilliard and Allen-Cahn equations [9], and it 
is the program we follow here for the MBE model. 

Our primary results are the following: for the equation 
of motion, Eq. ([T]), there exists a class of first order, semi- 
implicit steps 



ht+At = ht + At 



-V^ht - V ■ {{1 - \S/ht\^)Vht} 



biAtV^{ht+At - ht) + b2AtV\h 



t+At 



(2) 



that provides stable numerical integration for appropri- 
ate choice of the parameters bi and 62, as shown in Fig.[T] 
The results of our UvN stability analysis are presented 
as shaded regions while our numerical tests of gradient 
stability are plotted as points. Although UvN stability 
does not ensure gradient stability, we find that it is very 
effective in determining the gradient stable regions, both 
for single- and many-domain systems. The UvN stability 
conditions plotted here are independent of lattice type or 
details of the numerical method (e.g., finite difference ver- 
sus spectral methods). The difference in stability range 
for the single- and many-domain systems is revealed by 
the UvN stability analysis, which shows that the most 
unstable Fourier mode is that with its wavevector ori- 
ented with the local slope V/i. In the many-domain sys- 
tem, each mode samples many different slope directions, 
which acts to suppress the instability for the parameter 
range -1 < 61 < -1/2. 

Our results are consistent with Xu and Tang, who 
proved gradient stability for the parameters 61 < — 1 and 
62 = 1 Further work by one of us has led to a 

demonstration of gradient stability for the entire dark 
gray shaded region of Fig. [1] 17]. This proof will be pre- 
sented elsewhere, as it is considerably more general than 
the model considered here, and does not distinguish the 
single- and many-domain cases captured by the UvN sta- 
bility analysis. 

The remainder of the paper is as follows. In Sec. Ullwe 
review some of the properties of the model to provide nec- 
essary background for subsequent sections. In Sec. lIIII we 
present the UvN stability analysis, both for single- and 
many-domain systems. We describe the numerical tests 
of gradient stability in Sec. IIV[ as well as providing the 
details of our finite-different implementation of Eq. 
This is followed by a summary in Sec. [V] 



II. THE CONTINUOUS TIME MODEL 

In this section we provide motivation for the model we 
are considering and present some of its properties, show- 
ing in particular the instability to pyramid formation and 
the energy decreasing dynamics of the continuous time 
model. 




FIG. 1. (color online) Stability diagram for the parameters 
&i and &2 in Eq. ([2}. The UvN stable parameter values are 
shaded in gray, with the darker region corresponding to a 
single-domain system and the combined gray regions corre- 
sponding to a many-domain system. The points represent 
numerical tests of gradient stability: the (blue) circles are 
parameter values that are stable for single-domain systems; 
these together with the (purple) squares are stable for multi- 
domain systems; and the x are parameter values that were 
found to be unstable. 



The height field h{x,y,t), is defined in a co-moving 
frame so that its average is zero, and obeys a continu- 
ity equation. The current J has an equilibrium surface 
diffusion contribution equal to the gradient of the local 
curvature, Jsd ~ V(V^/i) and a non-equilibrium 
component Jne: 



dh 
'di 



— = - V • J = -V^h - V • Jt 



NE- 



(3) 



A noise term is omitted as this is considered to be irrel- 
evant for coarsening 8]. We consider the slope-selecting 
nonequilibrium current 



JNE = (i-|v;inv/i, 



(4) 



which gives Jne ^ V/i for small gradients, an uphill cur- 
rent due to the Ehrlich-Scliwoebel effect [3], and Jne = 
for slopes of unit magnitude. Inserting Eq. Q into the 
continuity equation ([3]) yields the equation of motion, 
Eq. dl]). 

Common variations on this model include slope- 
selecting currents that vanish for only a discrete set of 
V/i directions, reflecting the underlying crystalline struc- 
ture, and models without slope selection. The physical 
basis and experimental evidence for these various mod- 
els is described in [1, H, [T9l - [2l| and references therein. 
Material parameters have been absorbed into rescaling 
of lateral space dimensions, height, and time. 

The equation of motion, Eq. ([1]), can be written as a 
gradient flow 



dh 
'dt 



"6h 



(5) 
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for the free energy functional 



F[h] 



jd^x[i{vhf + \{i-m^f]. (6) 



Gradient flow results in a monotonically decreasing free 
energy, 



dt 



-F 



2 , 5F\ dh 



dh 

at 



<0. (7) 



As first noted by Eyre , the essential stability criterion 
for discrete time steps is to preserve the energy decreasing 
property of the continuous-time equation. This is known 
as the gradient stability condition. 

Next we review the the linear stability of the contin- 
uous time equation, which will be useful context for the 
von Neumann stability analysis in Sec. IIIII Consider a 
height field 



Hx,y,t) = + T]{x,y,t), 



(8) 



which consists of small deviations rj from a uniform slope. 
Inserting this into Eq. ([T]), linearizing in 77, and Fourier 
transforming to 7y(k, t) = J d^x exp(ik • x.)r]{x, y, t) gives 



dt 



= (fc2 - fc4 _ c^k^ - 2C^kl) ?7(k, t). 



(9) 



For an interface that is initially flat we set C = and 
obtain the growth rate for small fluctuations in the initial 
conditions: 



(10) 



step should be gradient stable, to preserve the energy- 
decreasing property of the continuous equation. How- 
ever, in this section we analyze instead von Neumann 
(vN) stability, i.e. the linear stability of the discrete step, 
Eq. ([2]). This analysis has certain advantages. It is rel- 
atively straightforward and, as shown in Fig. [T] and in 
Ref. [9^, it successful predicts the parameter range for 
gradient stability, as judged by numerical tests. Also, the 
method provides insight into the dynamics of the Fourier 
modes, which in the present case proves useful in clarify- 
ing the distinction between the single- and many-domain 
systems. 

We first present vN stability analysis on the Euler 
step, which results in conditional stability, i.e., a lattice- 
dependent upper bound on At. Then we consider our 
parametrized semi-implicit step and perform uncondi- 
tional vN stability analysis; that is, we seek parameter 
values which yield vN stable steps for any size At. Note 
that we will only impose vN stability on the equilibrium, 
sloped interface and not on the flat interface, where the 
linear instability is part of the physical behavior of the 
continuum equation. 

In addition to the time discretization, the spatial 
derivatives in our equation of motion must be treated 
by finite-difi^erence or spectral methods. Without spec- 
ifying the details of the scheme, we denote the Fourier 
transform of the two-dimensional numerical laplacian as 
A(k). In the continuum limit, A(k) — > — fc^. For spatially 
discretized systems, > A(k) > Amin, where the value of 
the lower bound Amin ~ — l/Ax^ depends on the details 
of the discretized laplacian. Our stability conditions will 
rely only on the universal upper bound of zero. 

We will use X{kx) to represent the Fourier transform 
of the numerical derivative second derivative d^/dx^. 



Long wavelength modes with fc < 1 are unstable and 
grow, which is exactly the instability that leads to pyra- 
mid formation. In the context of the Cahn-Hilliard equa- 
tion this is the spinodal instability Q. Note that the 
exponential growth of the mode is nevertheless accom- 
panied by a decreasing total free energy, as required by 
Eq. 0. ' 

For an equilibrium interface we set the slope C = 1 to 
obtain 



dyjKt) 
dt 



= -ik^ + kl)i^ik,t). 



(11) 



The negative right hand side indicates that height fluctu- 
ations about the equilibrium slope decay, and the uniform 
slope profile is stable. 



III. UNCONDITIONAL VON NEUMANN 
STABILITY ANALYSIS 

The goal in constructing a discrete time method is 
to be faithful to the physical behavior of the continu- 
ous time equation. In our case, this means our discrete 



A. Euler Step 

Our discrete time step, Eq. ([2]), reduces to an Euler 
step in the case 61 = 62 = 0. We plug in h ^ x + rj (i.e., 
slope C — 1), linearize in 77, and Fourier transform to 
obtain 



Vt+At 



1 + At 



{-A(k)2 + 2A(fc,)} 



(12) 



The vN stability condition is that the square bracket term 
has magnitude less than unity, to ensure fluctuations die 
away. The negative curly bracket term in Eq. (1121) has no 
lower bound in the continuum limit Ax 0, and thus 
the Euler step would be vN unstable for any size At. The 
situation is improved by the numerical derivative, which 
places a lower bound on the curly bracket terms, leading 
to vN stability for At < lAminl"^ ^ Ax'^. The analysis is 
essentially identical to what happens in the Cahn-Hilliard 
equation The Euler step provides an example of a 
lattice-dependent stability condition (relying on the lower 
bound of A(k) rather than the upper bound of zero) and 
it results in a fixed bound on the time step, regardless of 
the natural time scale of the dynamics. 
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B. UvN Stability for a Single Domain 

We return to our parametrized discrete step, Eq. ([2|), 
but now we leave bi and 62 unspecified. We seek to find 
ranges for the parameters wfiich will lift any restrictions 
on At, i.e., unconditional stability. We substitute Eq. ([8]) 
with slope C = 1 into Eq. ([2|), linearize, and Fourier 
transform. The resulting step can be written as 

[1 + Atcik)] 7i+At = [1 + A^7^(k)] fjt (13) 

with 

/:(k) = 6iA(k) +62A(k)2 (14) 

and 

7^(k) = 2X{k^) + biX{k) + (62 - l)A(k)2. (15) 

Before imposing the UvN stability, we note that it is 
necessary to have C(k) > so that the square bracket 
on the left of Eq. is non-vanishing for all At and k. 
This gives the requirement that 61 < and 62 > 0. 

Next, the UvN stability condition, |77t+At| < \Vt \ for all 
At and k, will be satisfied if C(k.) > |7?.(k)|. In the case 
that TZ(k) is positive, this gives the condition 

< £(k) ^7e(k) = -2A(fc^) + A(k)2, (16) 

which is intrinsically satisfied due to the non-positivity of 
X{kx). While here and below the k = mode saturates 
the bound, we can safely ignore it since it is static. 

The crucial condition, then, comes from imposing 
£(k) > —TZ(k), which becomes 

A(fc,) + 6iA(k) + (^62 - A(k)2 > 0. (17) 

The last term is positive for 62 > 1/2. Next, noting that 
A(fca,) > A(k), we have a lower bound on the remaining 
two terms: 

A(fc,) + 6iA(k) > (l + &i)A(k). (18) 

This will be positive provided that hi < —1. Thus, our 
conditions for UvN stability of a single-domain system 
are 

h < -1, &2 > 1/2, (19) 

which is plotted as the dark gray region of Fig. [1] 

Note that for bi slightly above —1, in the unstable re- 
gion, it is Fourier modes with X{kx) ~ A(k) that first vio- 
late Eq. (jl7|) . This means wavevectors k that are nearly 
oriented along the x-axis, i.e. the gradient direction of 
the equilibrium interface. 



C. UvN Stability for a Many-Domain System 

In a many-domain system, which is the relevant case 
for coarsening studies, we are not free to choose the coor- 
dinate axes to align the x axis with the interface gradient, 
since there are many facets with different gradient direc- 
tions. To analyze this case, we first about a single domain 
but with an arbitrary normal direction, parametrized by 
the polar coordinate 9 

h{x, y, t) = cos(6')x -I- sin(6')y + ri{x, y, t). (20) 

This follows through just as before, with the important 
stability condition Eq. p7)) becoming 

cos^eX{kx)+s\n^ eX{ky) + biX{\^)+[b2 - ^ A(k)2 > 0. 

(21) 

Now, if many domains are present in the system with 
essentially random orientations, then for any particular 
Fourier mode the above equation will be averaged over 
6*, giving {cos^B) = (sin^ 0) = 1/2. Using 

A(A:,) + A(fey) w A(k) (22) 

reduces Eq. ([2T|) to 

{b, + A(k) + {b2 - A(k)2 > 0. (23) 

Thus, our UvN stability condition for many-domain sys- 
tems is 

5i < -1/2, b2 > 1/2, (24) 

which is depicted as the combined shaded regions of 
Fig.[TJ The averaging over multiple orientations provides 
a greater parameter range of stability than the single- 
domain case. 

Note that in general Eq. ((22) is only an approximate 
relationship. It is a strict equality in the Ax — > contin- 
uum limit, and also in the common five-point stencil for 
the numerical laplacian on a square lattice, but for other 
choices of numerical derivatives it need not be exact. 



IV. NUMERICAL TESTS OF GRADIENT 
STABILITY 

Since the field equation of motion is nonlinear, von 
Neumann stability analysis is not sufficient to prove gra- 
dient stability. For that reason, we have conducted ex- 
tensive numerical tests for gradient stability for a range 
of 61 and 62 parameter values. We present the details of 
the numerical derivative implementation in an appendix, 
but we note here two important general features such 
an implementation should have. First, the local conser- 
vation law should be constructed to hold exactly, not 
just to some order in Ax, and second, the energy-decay 
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property of the continuous time equation should be main- 
tained when spatially discretizing. That is, the particular 
scheme of calculating the spatially discrete analog of the 
free energy F[h] in Eq. ^ and the equation of motion 
should be consistent, so that 



dF 



(25) 



is an exact relation, not just approximate to some order 
in Ax. 

For each bi and 62 value represented as a data point 
in Fig. [l]we performed the following tests. We evolved a 
512 X 512 sized lattices with lattice constant Ax = 1 out 
to a final time tmax- These systems were evolved using 
three different methods: an Euler step with At — 0.03 
out to a tniax = iC, a semi-implicit step with bi = 2.5 
and 62 = 1 and growing time step At = 0.03t^/^ out to 
time tmax = 10^, and the same semi-implicit parameters 
with a fixed time step At = 100 out to time i,„ax = 10®. 
For each of these cases we analyzed multiple runs and 
varied between random initial conditions and sinusoidal 
initial conditions with long and short wavelengths. 

At regular intervals during the evolution we tested a 
single step calculated via Eq. ([2]) with sizes varied be- 
tween 1 < A < 10^°. This step was used only for energy 
stability testing and did not contribute to the subsequent 
time evolution. Any time that the free energy was found 
to increase, that particular set of parameter values was 
identified as unstable. 

For the many-domain system, we used periodic bound- 
ary conditions and an initially flat interface (plus the ran- 
dom or sinusoidal fiuctuations). For the single-domain 
system, we first re-write the field equation of motion, 
Eq. ([T]) in terms of deviations from the uniform slope, 
giving 

^=-VS + 2dlv + 2d,\Vv\^ + 2{d.,r,)V^r, ^^^^ 

+ V-(|V77|2V7/), 

where dx = d/dx, and then constructed the analogous 
numerical implementation of this equation. This ap- 
proach was necessary to eliminate sensitivity to trunca- 
tion error. We imposed periodic boundary conditions on 
77, which corresponds to shifted periodic boundary con- 
dition on h. 

In Fig. [1] we show the results of this testing both for 
the single- and many-domain systems. The (blue) circles 
represent parameter values that were found to be stable 
for the single-domain system, that is, under all our test- 
ing, there were no single incidents of energy increase. The 
(purple) squares are parameters values that were found 
to be unstable in the single-domain system, but stable for 
the many-domain case. The remaining x are parameter 
values found to be unstable for both single- and many- 
domain systems. We find a striking degree of agreement 
between the predictions of UvN stability analysis and the 
numerical tests for unconditional gradient stability. This 
is one of our main results. 
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FIG. 2. (color online) The free energy density / — F\h]/L^y^ 
as a function of time, where the time evolution utilized a 
growing time step. At ~ t^^'^ . Simulation details are in the 
text. The expected t~^^^ decay is plotted for comparison. 



There is a small region for 62 < 1/2 where numerical 
tests find gradient stability. This can be understood from 
Eq. (|17l) as a lattice-dependent stability arising from the 
laplacian lower bound Amin- We have emphasized instead 
the lattice-independent stability boundaries, as these are 
more widely applicable. 

To illustrate the utility of these methods, we have sim- 
ulated the coarsening that results from an initially fiat in- 
terface, using a the stable step parameters 61 = —1.5 and 
&2 = I and a growing step size At — max(0.1, O.Oli^/'^). 
We performed 20 independent runs on a 2048 x 2048 lat- 
tice with Ax = 1, out to time tmax = 10*". 

Fig. [5] shows the decay of the free energy with time. 
Once equilibrated domains form, the free energy density 
F is proportional to the amount of edge in the system, 
which is inversely proportional to the characteristic size 
of the domains. Thus the free energy should decay as 
F ~ 1/L{t) ^ t^'^/'^. Our growing time step integration 
reproduces this result. 

Shown in Fig. [3] are snapshots of domain configurations 
for various times from a single run on a 512 x 512 lattice, 
with all other parameters as given above. 



V. SUMMARY 

We have parametrized a first order accurate discrete 
time step, Eq. ^ for MBE growth with slope selec- 
tion that, unlike the Euler step, is gradient stable for 
appropriate choices of the parameters bi and 62. We 
determined the stability range for these parameters via 
unconditional von Neumann stability analysis, and then 
tested these predictions with numerical tests for gradi- 
ent stability, as shown in Fig. [T] We find that the UvN 
stability analysis serves as an accurate proxy for uncon- 
ditional gradient stability, similar to the behavior of the 
Cahn-Hilliard equation [9|. 

Our stability analysis contained an implicit assumption 
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t = 1 000 



t = 5 000 





t = 25 000 



t = 100 000 



FIG. 3. (color online) Plotted is the laplacian of h{x, y, t), for 
a system evolved with a growing time step At ~ t^^^ . Simula- 
tion details are provided in the text. Positive values (troughs) 
are red, negative values (peaks) are blue, and the white re- 
gions are domains of uniform slope with zero laplacian. 



that the interface slopes do not exceed unit magnitude, 
which we justify by noting that the dynamics naturally 
select for this slope. This came into our UvN analysis 
by our choice to linearize about a unit slope domain. 
We note that the numerical tests for gradient stability 
contained no such assumption, so the agreement between 
the two approaches confirms validity of the unit slope 
assumption. 

Since UvN stability analysis is a relatively straightfor- 
ward process, we expect this approach to be generally 
applicable to a wide variety of phase field models. In 
particular, the case of slope-selecting MBE growth with 
preferred slope orientations, as studied by [5[ [l^, [1^, [13] 
and others, should follow identically to the case presented 
here. 
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Appendix: Finite Difference Scheme 

Here we present details of the spatial discretization 
scheme we used in our numerical tests. We present these 
in a discrete-space, continuous time picture, as our goal 
is to ensure that the conservative dynamics and the gra- 
dient flow are exact, i.e. preserved to all orders in Aa:. 
The essential condition for gradient flow is that the equa- 
tion of motion must be connected to a particular choice 
for the free energy functional such that 



dt 



d 



dh,. 



F 



(A.l) 



Local conservation is imposed by ensuring that the equa- 
tion of motion has the form 



dt 



1 r 
Ax 



{Jy}iJ+l/2 - {Jy}iJ-l/2 



(A.2) 



so that the same {Jx}i+i/2,j flows into /li+i.j and out of 
hi,j, and the same {Jy}ij+i/2 flows into /li.j+i and out 
of hi J. 

Our implementation uses an on-site finite-difference 
expression for V^/i, for which we take the standard five- 
point stencil, 

{V^/i}j,j = -^-^[hi+ij + hi-ij + hij+i + Hlj-i - 4/iij], 

(A.3) 

and the cell-centered expression for |V/ip, 



{iv/iin 



1+1/2J + 1/2 - 

+ {hi.j+i ~ hij)'^ + (/ii+i,j+i — hi+i,j)'^ 



(A.4) 



With these choices it is straightforward to show that 
d 



dh 



k,l 



k,l- 



(A.5) 



Our equation of motion is given by Eq. (jA.ip with the 
choice 



F 



(A.' 



By making use of Eq. (jA.Sp . the equation of motion can 
be shown to satisfy the discrete continuity equation (jA.2[) 
with current 

{Jx}i+l/2,j = {Jx^}i+l/2J + {Jx^}i+l/2,j (A. 7) 

where the surface diffusion current is 



i+l/2,J 



Ax 



(A.8) 
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and the nonequilibrium current is 



{■Jx h+l/2,j 



Aa 



1 - \{{\^h\^}^+l/2,J + l/2 + {|V/ip},+ i/2,,_l/2) 

(A.9) " 



and analogous expressions for {Jy}ij+i/2- The discrete 
from of the free energy, Eq. (jA.6l) . was used in aU our 
numerical tests for gradient stability. 
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